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Abstract 

Experiments suggest that localization via self-trapping plays a central role in the behavior of 
equilibrated low mass particles in both liquids and in supercritical fluids. In the latter case, the 
behavior is dominated by the liquid-vapor critical point which is difficult to probe, both experimen- 
tally and theoretically. Here, for the first time, we present the results of path-integral computations 
of the characteristics of a self-trapped particle at the critical point of a Lennard-Jones fluid for a 
positive particle-atom scattering length. We investigate the influence of the range of the particle- 
atom interaction on trapping properties, and the pick-off decay rate for the case where the particle 
is ortho-positronium. 
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The system consisting of a massive particle equilibrated in a host fluid, known to physicists 
and chemists as Brownian Motion, has played a seminal role in the development of Statistical 
Physics 0, ^|. However, the opposite regime of an equilibrated low mass particle is equally 
challenging and manifests a richer set of behaviors [3|. Except at very high temperatures, 
quantum mechanics is required to model the low mass particle. Possible quantum particle 
(qp) candidates are an electron, positron, or positronium atom, while the host can either 
be a dense gas or liquid below the critical temperature, or a supercritical fluid above it. 
Experimental measurements of the properties of a low mass particle equilibrated in a fluid 
strongly suggest that it can induce a local, mesoscopic, deformation in the fluid in which the 
qp becomes self-trapped, or localized. ^, ^ | Since the qp has a long deBroglie wavelength, 
intuitively we anticipate that it simultaneously interacts with a large group of atoms or 
molecules in the host fluid forming a mesoscopic region of altered fluid density. Depending 
on whether the effective qp-atom interaction is attractive or repulsive, the local density of 
the host is either augmented or suppressed near the qp, resulting in the formation of either a 
"microdroplet" or " microbubble" . The intuitive picture is completed by imagining that the 
qp occupies the ground state of the potential well induced by the formation of the density 
inhomogeneity, i.e. the droplet or bubble, thus stabilizing the deformation. 

The positron-atom interaction is characterized by a negative scattering length so the 
experimental manifestation of self-trapping is a decrease in the positron lifetime due to 
the increase in the local electron density resulting from the formation of a droplet. The 
reverse is true for ortho-positronium. Angular momentum conservation eliminates the two 
photon decay process for the triplet state. As a result of the long natural lifetime of o-Ps 
in the vacuum, about 142 ns [4], its positron can annihilate more readily with an atomic 
electron of the host. A consequence of the fermionic repulsion between the electron in o-Ps 
and the electrons of the host fluid is that this "pick-off" annihilation rate is reduced when 
self-trapping occurs, demonstrating bubble formation [4]. Depending on the choice of the 
host, the electron- atom scattering length can have either sign, so each behavior is possible. 
While the most dramatic manifestation of self-trapping of the positron or positronium is 
a significant deviation from linearity in Arrhenius plots of the decay rate versus average 
density, the signature of electron self-trapping is a change in mobility [5|. 

In addition to the liquid state jf|, self-trapping occurs in a broad, super-critical region 
of density, p, and temperature, T, surrounding the liquid-vapor critical point (p c , T c ) 



region 

m 
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Since the isothermal compressibility diverges at the critical point, this is not surprising: 
The qp can more easily alter the local density in this region of temperature and density. 
However, as a result of the large density gradient induced by the earth's gravitational field, 
there are few reliable experimental studies of self-trapping close to the critical point 0. 
Thus, although it dominates the self-trapping regime above TV, the effect of close proximity 
to the critical point on self-trapping is not generally known 8|. 

The theory of self-trapping has evolved through different stages: In the earliest mod- 
els, the qp occupies the ground state of a spherical step potential which is assumed to be 
proportional to the local fluid density Q]. In modelling the qp-host interaction the atomic 
nature of the host is ignored, and it is simply represented by a type of jellium. By mini- 
mizing the free energy of the cm-fluid system, it is poss.b.e to show that the deformation 
is stable in a bounded region of temperature and density near the critical point j9(. Later 
versions of mean field theory (MFT) permitted a continuous density profile ]?], Q| and took 
into account the atomic nature of the host at an intermediate level An interesting 

improvement was obtained by employing the Percus-Yevick equation to include the effect 
of qp-atom and atom-atom correlations in the mean field formalism [3]. An advantage of 
MFT is that computations are reduced to numerically integrating a pair of coupled ordi- 
nary differential equations. However, in practice, they have only proved useful for fluids at 
low temperatures 

0,0, HQ- 

This problem may arise because mean field theories only 



include a single bound state for the qp. More microscopically complete models have evolved 
during the last few decades based on the Feynman/Kac path integral ^| which overcome 
the major shortcomings of the earlier work. These account for the details of the qp-atom 
interaction and implicitly take into account local fluctuations in the disturbed fluid and the 
state of the qp. Two approaches have been em ploy ed, one based on an approximate analytic 
model requiring an educated choice of closure Ljjllal) an d the more direct alternative using 

nn 

Monte Carlo algorithms [13J, Hj| . While the latter avoids approximation, the computational 
cost is greater. 

Here we employ path integral Monte Carlo (PIMC) to investigate self-trapping at the 
liquid- vapor critical point. We take advantage of recent improvements in fluid equilib- 
rium theory which give accurate critical point parameters for the Lennard- Jones (6,12) 



potential. 



18l ? ] We model the qp-atom interaction with a hard-sphere potential and study 



the dependence of the physical and statistical properties of self-trapping on the hard sphere 
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diameter, Rh s . In particular, we investigate the dependence on Rh s of both the spreading of 
the qp wavefunction and its influence on the local deformation of the fluid. Since a repul- 
sive qp-atom interaction is representative of positronium [dj], we also estimate the pick-off 
annihilation rate for a simplified model of the atomic charge distribution. 

Although no experimental measurements of the pick-off decay rate have been carried out 
precisely at the critical point, for the case of Xenon they exist for two supercritical temper- 
atures, 300-K" and 340^, over a large density range and strongly suggest the existence of the 
Ps self-trapped state [13 . In previous work we have used PIMC to carry out simulations 
of Ps in Xenon at these temperatures, so it is natural to select the critical point of Xenon as 
a test case and a mass of 2m e for the qp. The truncated version of the Lennard- Jones 6-12 
potential was chosen to represent the inter-atomic potential. Wilding has established numeri- 
cally accurate connections between the L- J distance and energy parameters (a and e) and the 
critical temperature and density (T c and p c ), namely p c a 3 = 0.3197 and kT c /e = 1.1876 [18]. 
The experimental values of the critical density and temperature of Xenon are T c = 289i^ 
and p c = 5.299 x 10 -3 atoms/A 3 yielding a = 3.92A and e/k = 243.5K. 

The computational simulation of the qp-fluid system using a complete quantum me- 
chanical description would be intractable; thus approximate representations of the system 
are required. Except at very low temperatures, the translational degrees of freedom of an 
atomic fluid can be approximated with classical mechanics. Thus the qp-fluid system is well 

represented by a hybrid classical-quantum Hamiltonian: 

n fc2 *r 



F = ^Pj/2M + E tf(|H*-H*l) + ff,p, ff9P = -2^A + XV(|r-R,|). W 



Here Pj and R., are the momenta and positions of the N host atoms, r is the position of the 
qp, while U and V are the pairwise additive atom-atom and qp-atom interaction potentials. 
This is known as the adiabatic formulation and results in a hybrid partition function, (^(J 



where here R represents {Rj}, the complete 3iV dimensional configuration space of the 
atoms with classical potential energy U(R), and A is the atomic thermal wavelength. Thus, 
in principle, the quantum statistical average of the physical operator O can be computed 
from 



N>j>k>l 
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To obtain a formulation which is useful for computation, we follow the Feynman-Kac path 
integral construction First, by applying the Trotter formula to Z, we can express the 
trace over r as a sum over discretized paths of P steps. For sufficiently large P the kinetic 
and potential energy operators approximately commute in each step, yielding the following 
expression for the partition function; 

Z P =f dRe-^n / drie-^-^^piTuT^B.P/P)), (4) 

o 

where A is the qp thermal wavelength. If we take p to be exp(— j3V/ P), equation (4) is known 
as the primitive approximation. Thus, in the discretized Feynman-Kac path integral jl^j ]. 
the qp is represented by a closed chain of P classical "pseudo-particles", or "slices" in imagi- 
nary time ti = (i/P)/3h, with harmonically coupled nearest neighbors. Each pseudo-particle 
in the chain interacts with each atom through the potential V/P.The chain is equilibrated 
in a fluid at the augmented temperature PT jl?! I^ . In the limit P — > oo the correspon- 
dence is exact. The spread of the chain corresponds to the uncertainty in position of the qp. 
The equivalence allows Monte Carlo methods developed for classical systems to be used to 
compute quantum mechanical equilibrium properties. Path integral Monte Carlo has been 
used successfully by many groups, including our own, to compute the equilibrium properties 
of quantum systems J3]. 

To insure convergence, while maintaining a constant temperature and average atomic 
density, the system size, the number of chain-particles P, and the number of statistical 
samples are increased until there are no significant changes in the calculated equilibrium 
properties. The algorithm we employed here depends on five parameters, the density, tem- 
perature, Rhs, the number of fluid atoms, N, and P. Convergence was improved by using 
an image potential to smooth out the singularity of the hard sphere, qp-atom interaction 
20, Q- A complete description of this method is given in our previous publications (see 
20 1 and references cited within). For this initial study we selected three values of Rh s ; 0.5A 



, 5.0A and 9.5A while P was chosen to be 2000 as in our previous computations 20|. Since 
most of the cpu time was expended on repositioning the atoms, the large value of P did not 
pose a problem and insured convergence. To minimize edge effects, the length of the box in 
which the simulations were carried out was set to three times the qp thermal wavelength. 

Important characteristics of the qp are the root mean square displacement (D(t)) between 
two particles along the chain separated by imaginary time t, and the density of chain particles 



a distance r from its center of mass (pcm(r)). D and p cm are both direct measures of self- 
trapping. In an extended state the displacement between particles separated by t increases 
with t, reaching a maximum at t = -P/2, according to a unimodal parabolic distribution. 
However, in a self-trapped state the chain is highly confined and thus, except for values of 
t near the end points at and P, D(t) is roughly constant. The function p cm (r) represents 
the mean probability density of the qp wavefunction, and becomes more peaked around 
the chain com as the chain becomes confined within its self-trapped bubble. A detailed 
discussion of these quantities and how they are computed can be found in our previous 
paper on self-trapping of Ps in Xenon j^jj • 

Figured] is a plot of D versus t for the three Rh s values considered here. As Rh s increases, 
the figure shows that the shape of the plots changes from parabolic, corresponding to a 
nearly free particle, to one that is essentially constant in the central region, indicating that 



the chain becomes more confined 
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2l|. Increasing the value of Rh s expels fluid atoms 



from the vicinity of the chain particles into regions which don't overlap with them. At the 
same time the pressure applied by the fluid atoms compresses the chain. Since there are 
no repulsive forces between the pseudo-particles comprising the chain, their density near 
the chain com can become quite large. This can be seen in Figure 121 a plot of p cm versus 
distance from the com, for the same Rh s values shown in Figure ^ 

Information concerning the local deformation of the fluid induced by the qp is provided by 
the chain-fluid and com radial distribution functions, gj p (r), g/ cm (r), yielding, respectively, 
the mean local density of fluid atoms a distance r from a chain pseudo-particle, and the 
chain center of mass. As usual, they are normalized to unity in the large r limit. In earlier 
work we showed that the former can be used to directly determine the o-Ps pick-off decay 
rate |2Q], while the latter provides direct evidence of self-trapping. In an extended state 
fluid atoms are able to penetrate into the vicinity of the chain com indicating that the qp is 
relatively spread out. On the other hand, in a completely self-trapped state, ground state 
dominance prevails, the chain is folded upon itself within the volume of the bubble, and 
fluid atoms are totally expelled from this region. gf cm {^) can be used to directly compute 
the average number of fluid atoms excluded from the vicinity of the chain. Figure El a plot 
of gfcm versus position, indicates the strong exclusion of the fluid atoms from the vicinity 
of the chain at the largest value of Rh s - The significant value of p cm at the origin for the 
smallest value of Rh s provides evidence that the atoms can penetrate the trapping region 
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and suggests that the qp is in an extended state. This is supported by the parabolic shape 
of D(t). Hints of slow oscillation on large scales, which we are unable to resolve with the 
present system size, can be seen in Figure El An open question is whether these are simply 
finite size effects or represent a coupling between the correlation length of the fluid and the 
thermal wavelength of the qp. In our earlier work jscj short range oscillations occurred in 
Pcm over distances on the order of Rh s at p = 2p c , indicating stacking of the fluid atoms in 
layers around the qp com. We see from the figure that at the critical point they are absent. 

From g f cm we found that the volume of the excluded fluid atoms from the trapping region 
is at least twice as great for Rh s = 9.5A than for either 0.5A or 5. OA. Since increases in 
Rh s will lead to decreases in the number of fluid valence electrons near the qp available for 
annihilation, for the case where the qp is o-Ps the pick-off decay rate is expected to decrease 
with increasing Rh s - We computed the decay rate at the three values of Rh s considered 
above (see table QJ). Taking ratios, we found an approximately exponential decrease with 
Rhs with a characteristic length of 1.61 A . In this model the Ps atom is treated as a 
composite particle in its internal hydrogenic ground state with an exponentially decreasing 
wavefunction representing the spread of the positron from the Ps center of mass. This 
is justified by the small ratio, k 4 x 1CT 3 , of kT c to the positronium excitation energy. 
For simplicity, the electron density around the host fluid atoms is modelled as a delta 
function. Thus the decay rate is dependent upon both the portion of the positron ground 
state wavefunction that manages to leak beyond Rh s and the number of fluid atoms centers 
it overlaps. The exponential decrease in the decay rate with increasing Rh s results from the 
decrease in both the density of fluid atoms and the amplitude of the positron wavefunction 
in the shell surrounding the hard sphere surface. 

In contrast with mean field models, the power of PIMC is that it reveals the complete 
picture of both the distribution over qp quantum states and the response of the fluid to the 
qp " impurity . In genera, we havefcund that eompared with a liquid [6] and our own PIMC 
computations at higher density |2QJ], the mesoscopic region in which the qp is localized is 
much larger at the critical point. Moreover, the density profile of the fluid has a different 
shape - the "walls" are not nearly as steep and the structure is less sharply defined, i.e. it is 
soft. Our PIMC calculations for a supercritical fluid with p = p c , but at hi ghe r temperature, 
show similarities with the critical point, but the behavior is less extreme (20 1. 

In future work we plan to employ PIMC to directly evaluate the density of states. This 
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TABLE I: Natural log of the decay rate for three hard sphere diameters. 



Rhs (A) Ln(Decay Rate) 
0.5 -3.0 
5.0 -20.0 
9.5 -37.0 



will enable us to determine the angular distribution of the annihilation photons from para- 
positronium, which can be directly compared with experiment. To date this has only been 



a. 



It will be interesting 



approached with a semi-empirical formulation of mean field theory 
to see what changes result from an ab initio computation which completely includes the 
effects of correlated fluctuations. 
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FIG. 2: Density of chain particles as a function of distance from the chain centroid, p cm (r), for the 
smallest (0.5A), average (5.0A), and largest (9.5A) values of Rh s - 
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FIG. 3: Density of fluid atoms from the chain center of mass, gfcmif), versus position for the 
smallest (0.5 A), average (5.0 A), and largest (9.5 A) values of Rh s - 



12 



